Dynamics of Saturated Energy Condensation in Two-Dimensional Turbulence 
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In two-dimensional forced Navier-Stokes turbulence, energy cascades to the largest scales in the 
system to form a pair of coherent vortices known as the Bose condensate. We show, both numerically 
and analytically, that the energy condensation saturates and the system reaches a statistically 
stationary state. The time scale of saturation is inversely proportional to the viscosity and the 
saturation energy level is determined by both the viscosity and the force. We further show that, 
without sufficient resolution to resolve the small-scale enstrophy spectrum, numerical simulations 
can give a spurious result for the saturation energy level. We also find that the movement of the 
condensate is similar to the motion of an inertial particle with an effective drag force. Furthermore, 
we show that the profile of the saturated coherent vortices can be described by a Gaussian core with 
exponential wings. 

PACS numbers: 47.27.-i, 47.27.De, 47.27.E- 
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I. INTRODUCTION 

Two-dimensional (2D) hydrodynamic turbulence is 
fundamentally different from its three-dimensional coun- 
terpart. In 2D, small vortices can merge to form bigger 
coherent vortices. This is because the equations of ideal 
hydrodynamics in two dimensions have, in addition to en- 
ergy, also enstrophy as a conserved quantity. With an ex- 
ternal force at an intermediate scale and viscous dissipa- 
tion, energy inversely cascades to larger length scales and 
enstrophy directly cascades to smaller length scales [1-3] . 

Let us now consider 2D turbulence in a finite domain of 
size L. The smallest wave number allowed in this system 
is ki = 2iz/L. Due to the inverse cascade, energy piles 
up at k± provided there is no large-scale friction. This 
phenomenon, sometimes called Bose condensation in 2D 
turbulence (see Fig. 1), was first predicted by Kraichnan 
[1] . It was studied numerically by Hossain et al. [4] , Smith 
and Yakhot [5, 6], Chertkov et al. [7] and experimentally 
by Paret and Tabeling [8] , Xia et al. [9] . 

Following standard convention, we refer to the modes 
at |fe| = ki as the condensate in this paper. For a fixed 
non-zero viscosity, v, the energy of the condensate vor- 
tices cannot grow without limit but saturate [10, 11]. The 
saturation occurs at time scales of the order of \jvk\. 
This is an unusual example of viscous effects playing an 
important role in fluid turbulence at large length scales. 
In this paper we show, by direct numerical simulations 
(DNS), that the saturation of the condensate energy is 
not only determined by viscosity but also by the forcing 
wave number k\. Furthermore, we will demonstrate that 
the direct enstrophy cascade must be well resolved for 
accurate numerical determination of the saturation. 



Motivated by the analogy between the formation of 
large-scale structures in two-dimensional turbulence and 
the large-scale dynamo process in three-dimensional he- 
lical magnetohydrodynamic turbulence [12], we propose 
a simple three-scale model which is able to capture the 
important aspects of our numerical results. We further 
show that the Lagrangian dynamics of the condensate 
vortices can be described by Langevin equations for par- 
ticles with inertia. Finally, we measure the profile of the 
saturated coherent vortices, which consist of the conden- 
sate and its higher harmonics. The vorticity at the cores 
are well fitted by a sharp Gaussian; while the wings fall 
off exponentially. 
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FIG. 1: (Color online) Pseudocolor plot of vorticity showing 
Bose condensation. Red (the top-left vortex) and blue (the 
bottom-right vortex) represent positive and negative vorticity 
in physical space, respectively. The color scale is shown in the 
color bar on the right, which is chosen to make the fluctuation 
visible. The vorticity is normalized so that max |cj = 1. 
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al6 


1CT 3 


16 


80 


500 


405 


512 2 


483 


323 


0.80 




a24 


10~ 3 


24 


91 


500 


356 


512 2 


520 


281 


0.79 




a32 


1CT 3 


32 


100 


500 


308 


512 2 


582 


234 


0.76 




bl6 f 


5 x 10~ 4 


16 


113 


1000 


891 


512 2 








single- precision, crashed 


bl6 


5 x 1CT 4 


16 


113 


1000 


891 


512 2 


922 


696 


0.78 




b24 


5 x 10~ 4 


24 


129 


1000 


832 


512 2 


947 


643 


0.77 




b32 


5 x 1CT 4 


32 


142 


1000 


772 


512 2 


970 


580 


0.75 




cl6 


2 x 1CT 4 


16 


178 


2500 


2373 


1024 2 


2223 


1818 


0.76 




c32 


2 x 10~ 4 


32 


224 


2500 


2226 


1024 2 


2287 


1714 


0.77 




c64 


2 x 1CT 4 


64 


283 


2500 


1928 


1024 2 


2409 


1385 


0.72 




c64 b 


2 x 10 -4 


64 


283 


2500 


1928 


512 2 


1985 


791 


0.41 


converge to wrong answer 


dl6 


10" 4 


16 


252 


5000 


4859 


1024 2 


4145 


3581 


0.74 




d32 


10" 4 


32 


317 


5000 


4691 


1024 2 


4006 


3349 


0.72 




path 1 


5 x 10~ 4 


32 


142 


1000 


308 


512 2 




580 


0.75 


restart from saturated b32 



TABLE I: List of simulations: The first 13 simulations are used to verify our three-scale model (16). The last row describes 16 
restart-runs which are initially in saturated states. They are used to verify the properties of Ar . The details of the simulations 
are described in Sec. II. 



II. NUMERICAL SIMULATIONS 

Let ip be the 2D (scalar) stream function. The velocity 
is then u — {d y ip, —d x tp) and the z component of the 
vorticity is co = — V 2 V>- We solve the 2D incompressible 
Navier-Stokes equations in the vorticity-stream-function 
formulation, that is, 

d t <jJ — J(V>,w) = vV 2 lj + g, (1) 

where the Jacobian determinant is given by 

J(1>,w) = (d x ip)(d y u) - (d x Lo){d y ^). (2) 

We use periodic boundary conditions with L — 2tt so 
fci = 1, and denote the Fourier transform of to by cj^- 

In Eq. (1), g is the z component of the curl of an ex- 
ternal force. Using ( • ) to denote ensemble averages, the 
Fourier transform of the force, gk, is taken to be random 
and white-in-time with zero means and variance 

{gl{s)-g k {t))=f 2 k 2 5{t-s). (3) 

It is achieved by choosing a random phase for g^ at each 
time step. In the above expression, f[ is the forcing am- 
plitude and ki is the forcing wave number. To ensure 
isotropy, we select k randomly in a shell with radius fei 
and then round it off to the nearest grid point in Fourier 
space. The effective width of the shell is approximately 
k\. For this force, the average rates of energy and en- 
strophy inputs are respectively f 2 and f 2 k 2 . 

We use a spectral Galerkin scheme in space and a low- 
storage third-order Runge-Kutta/Crank-Nicolson time 
stepping scheme [13]. The time steps At are chosen with 
a Courant number of 0.5. The non-linear term is treated 
explicitly and the diffusion term is treated implicitly. The 



stochastic forcing is integrated by the Euler-Maruyama 
method [14]. We use NxN grid points with the Galerkin 
cutoff at k G = [(N - 1)/3J + 0.99, where |_ ■ J denotes the 
floor function. The value 0.99 is chosen such that the 
comparison k 2 < fc^ is accurate enough even in single 
precision with N <~ 2048, although almost all simulations 
in this paper are done with double precision. 

Our code is implemented in CUDA C and runs on graph- 
ics processing units (GPUs), which are massively paral- 
lel "stream" processors. With an nVidia Tesla C2050 
graphic card, our code is over an order of magnitude 
faster than codes running on a single CPU core. Be- 
cause of communication overhead, our code outperforms 
distributed-memory parallel codes running on 32 cores. 
This speed-up allows us to integrate the problem over 
very long code time and study the saturation of the con- 
densate vortices. We have released our Spectral Galerkin 
2D code (sg2) as an open-source project under the GNU 
General Public License version 3. The project is hosted 
by the Google Code http : //sg2 . googlecode . com/. 

III. RESULTS 

We have performed a series of simulations as shown in 
Table I with different sets of input parameters v and k\. 
In all of them, the initial conditions are Uk(t = 0) = 0. 
We set the forcing amplitude at /; = 1 and choose the 
forcing wave number fei so that k\ <C fei <C fed- The 
dissipation wave number fed is given by 

fed = ^-1/2, (4) 

where r) is the forward enstrophy flux. Further details of 
the runs are given in Table I. 
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FIG. 2: (Color online) Unnormalized energy evolution in sim- 
ulations al6 - d32. The thick red (lowermost), green, blue, 
and black (uppermost) curves are for different viscosity v (see 
labels); while solid, dashed, and dotted styles denote different 
forcing wave number ki (see legend). The green diamonds and 
blue squares are for simulations bl6 f and c64 b , respectively. 
They correspond to single precision and under-resolved sim- 
ulations. 



In our simulations, energy inverse cascades to small 
Fourier modes till fej and forms a condensate at this 
mode. We plot the total energy E{t) for simulations 
al6 - d32 in Fig. 2. The different colors and line styles 
represent different choices of viscosity v (see labels) and 
forcing scale k\ (see legend). The blue squares and green 
diamonds correspond to a single precision (bl6 f ) and an 
under-resolved (c64 b ) simulations, which we will com- 
ment in Sec. IV. For intermediate time, the growth of 
energy of this condensate is consistent with earlier re- 
sults [7]. Nevertheless, the asymptotic E(t) becomes in- 
dependent of time in all the cases. This saturation hap- 
pens at a time scale determined by the viscous time at 
the largest length scale (i.e., t„ = 1/vkf). This requires 
very long integration times and hence has not been ex- 
plored by earlier simulations. We are able to reach such 
late times by virtue of using a GPU code. 

At later times, most of the total energy comes from 
the energy of the condensate. Further understanding of 
the growth and saturation of the condensate can be ob- 
tained by studying the energy spectrum. We compute the 
shell-integrated energy spectrum Ek in simulation c32 at 
t = 1000 using a fixed bin size Ak = k\. The result is 
plotted as the thick solid curve in Fig. 3. The condensate 
is well developed but not fully saturated. The dashed and 
dotted lines in the figure are proportional to /c -5 / 3 and 
fc~ 3 , respectively. Although the inertial ranges are nar- 
row, the inverse energy cascade and forward enstrophy 
cascade are consistent with spectral slope —5/3 and —3. 
However, the condensate strongly deviates the inverse 
cascade spectrum at large scale. This additional feature 
motivates us to introduce a simple three-scale model. 




1 10 100 

Wavenumber k 

FIG. 3: Energy spectrum in simulation c32. The thick solid 
curve is the energy spectrum Ek computed at t — 1000 with 
bin size Ak — ki. The dashed line shows k~ 5 ^ 3 ; while the 
dotted line is proportional to k~ 3 . The gray areas, which 
represent our three-scale model, are drawn in proper scales. 
The condensate at zone (i) contains a large amount of energy. 
The broken power law in zones (ii) and (iii) have slopes —5/3 
and —3, respectively. See Sec. Ill A for details of the model. 



A. Energy and enstrophy evolution 

Given the energy spectrum Ek, the total energy, en- 
strophy, and palinstrophy can be computed by the one- 
dimensional integrals E = J E^dk, Z = J k 2 Ekdk, and 
P = J k 4 Ekdk, respectively. Together with the average 
energy and enstrophy input, the evolution of these quan- 
tities is governed by the following equations 

d t E = -2vZ + ft, (5) 
d t Z = -2vP + tfkf. (6) 

Our three-scale model is constructed in the following 
fashion. We refer to the corresponding wave number 
ranges as zones (i) through (iii) as shown in Fig. 3. (i) At 
the wave number of the box k\, we have the conden- 
sate modes, which contain most of the energy at satu- 
ration. We use the shorthand E(t) to denote the time- 
dependent energy at this Fourier mode, (ii) In the range 
k\ < k < ki, the energy inverse cascades from the for- 
cing wave number ki and forms an Ej~ ~ fc~ 5 / 3 spectrum. 
Note that the enstrophy spectrum Zk — k 2 Ej~ peaks at 
the forcing wave number fcj. (iii) The forward enstrophy 
cascade in the small scales gives the spectrum Ek ~ fc~ 3 
with a sharp cutoff at |fe| = fe<j- It is necessary to as- 
sume this cutoff to prevent the ultraviolet divergence of 
the total enstrophy and palinstrophy, although not total 
energy. 

Let us now denote by E the energy in the \k\ = k\ 
mode (i.e., the condensate) and by E 1 the energy of the 
rest of the system (i.e., E = E + E'). Following standard 
mean-field theory, we call them, respectively, the mean 
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and fluctuating parts, even though E' contributes to the 
coherent structure as we will show in Sec. IIIC. Similar 
notations are applied to the enstrophy and palinstrophy 
too. We can then split the energy equation (5) into two 
parts, namely, 



d t E = -2vk\E + e, 
d t E' = -2vZ' + ft -e, 



(7) 
(8) 



where e is the inverse energy transfer rate to the conden- 
sate. Similarly, with the help of Eq. (7), we can rewrite 
Eq. (6) as 



d t Z' = -2uP' + tfkf - kfe. 



(9) 



The last term kfe describes the enstrophy transfer rate to 
the condensate. Its existence implies that, although most 
of the enstrophy is transferred toward smaller scales, 
there is a small leakage towards the largest scale. 

Equations (7) through (9) are exact but not closed. 
Hence, it is not yet possible to solve them simultane- 
ously. We parametrize the fluctuating palinstrophy and 
enstrophy by 7 and T such that 



P' = P- k\E = yk%Z', 
Z' = Z- kfE = TkfE'. 



(10) 
(11) 



Substituting the above equations into Eqs. (8) and (9) we 
obtain two simultaneous equations. We then eliminate e 
between them to obtain a dynamical equation for Z', 
which can be solved to obtain 



Z'(t) = Z' x 



1 



-{t-t )/T z , 



(12) 



where we have used the initial condition Z'(to) — and 



Z' 



k\ 



2v -tkl 



with characteristic time scale 



1 rfc? 



k\ 



2vTk? 7 fc2 - k\ ' 



(13) 



(14) 



Now note that the enstrophy saturation time scale 
Tz' ~ l/2^7fcj is substantially faster than the energy sat- 
uration time scale, which is of the order of vk\. Hence, 
for simplicity, we can replace e by its value at late time 
in Eq. (7), which can be obtained by replacing Z 1 by Z^ 
in Eq. (8), 



eoo = lim e(i) = /; 

t— >oo 



ML = fi 



kf 



(15) 



Substituting this back into Eq. (7) and solving for E, we 
obtain 



E{t) = E c 



1 



,-(t-*o)A 



(16) 
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FIG. 4: (Color online) Normalized energy evolution in simu- 
lations al6 - d32. The thick red, green, blue, and black curves 
(all under the "Simulations" label) are for different viscosity 
(see also labels in Fig. 2); while solid, dashed, and dotted 
styles denote different forcing wave number fci (see legend). 
We normalize the energy and time using the results from our 
three-scale model. Almost all the curves collapse onto each 
other. The gray solid curve, green diamonds, and blue squares 
are the theoretical prediction, b!6 f , and c64 b , respectively. 



where we have used the initial condition E(to) = 0. This 
solution saturates to 



- f 2 



7*2-*? 



2vk\ 7*^ - 
with a characteristic time scale 

Tg = l/2vk\. 



k\ 



(17) 



(18) 



In the saturation level (17), the first part J? /2vk\ can 
be easily derived by dimensional arguments; while the 
correction factor (7*^ — kf)/(-fk"^ — k 2 ) < 1 comes from 
our three-scale model. This result is consistent with the 
energy bounds derivation by Eyink [10]. 

To obtain the values of Eoo> we need to solve for *d 
and 7. For the former, we combine definition (4) and the 
forward enstrophy flux in Eq. (9) to obtain 



f- /3 (kf-kf) 



2\l/6, -1/2 



(19) 



For the latter, we assume that in the inverse cascade 
Ek ~ fc -5 / 3 , and in the direct cascade Ek ~ fc -3 . These 
two power laws must match at fe = fej. Using these, we 
can compute the dimensionless number 



7 



3 + 41n(* d /fci) 



< 1. 



(20) 



With this, the model no longer has any free parame- 
ter. In Fig. 4, we re-plot E(t) with the vertical axis nor- 
malized by Eqo and also scale the time axis by 7V. This 
collapses all the different time traces, thus providing sup- 
port to our three-scale model. Note that the non-trivial 
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correction factor given in Eq. (17) is necessary to obtain 
this collapse. Note also that, the expressions for Eoo and 
Z'^ (but not the time scales) can be derived directly by 
demanding that we reach a stationary state at late times. 
Imposing the conditions d t E = d t Z = to Eqs. (5) and 
(6), Eqs. (17) and (13) follow immediately. 

Although all the results from different simulations col- 
lapse onto each other, this collapsed curve lies systemat- 
ically below the analytical solution shown as light gray 
curve in Fig. 4. To estimate this discrepancy, we fit the 
numerical curves by the form of Eq. (16) to find numer- 
ical values of the total energy and the saturation time 
scale. These values are compared with their theoretical 
prediction in Table I. The time scales agree to a good 
accuracy. The ratios between the theoretical Eoo and the 
fitted Eoo are listed in the tenth column of Table I. The 
average of these ratios (not including bl6 f , c64 b , and 
path 1 ) is ( = (Eoo/ Eoo) ~ 0.76. We speculate that this 
systematic difference exists because we have ignored the 
feedback of the condensate to the spectrum in our three- 
scale model. We will discuss this speculation with more 
details in Sec. IV. 

Comments on the remaining run in Table I are now 
in order. The blue squares are for the c64 b run. Al- 
though they look perfectly reasonable, the resolution of 
the simulation, 512 2 , which gives fcc ~ 170, is too low 
to resolve the actual enstrophy dissipation wave number 
k d = 283. Because spectral Galerkin methods have very 
little numerical dissipation, the forward cascaded enstro- 
phy cannot be removed, which artificially increases Z'^ 
and decrease eoo- Hence, the steady state energy Eoo 
converges to an (incorrect) lower value. In the language 
of numerical analysis, this is an inconsistent numerical 
solution. 

The green diamonds in Fig. 4 are for the bl6 f run, 
which uses single-precision floating-point numbers to rep- 
resent Wfc instead of double precision. The numerical so- 
lution behaves properly at the beginning, but starts to 
diverge from the true solution around t — 1000 = t^. 
This is not surprising because the relative round-off error 
of single-precision numbers is about ~ 10~ 7 . A random 
walk of round-off error in Wfc leads to a linear grow of er- 
ror in the energy. Indeed, there are roughly 10 7 steps by 
the time we reach t = r-g, which gives the order of unit 
error seen in the figure. To date, a significant fraction of 
scientific computing in CPUs are done with single preci- 
sion. This turns out to be insufficient for our purpose. 



B. Movement of Condensate Vortices 

After studying the saturation level and time scale, we 
focus on the movement of the condensate vortices at late 
times. The last row in Table I represents 16 different sim- 
ulations path 85 - path 100 . The superscripts in the names 
have the following meaning. For path 1 , we pick the i-th 
output file of b32 and restart the simulation each time 
with a different realization of the random force. In other 



words, for the i-th run the random number generators 
are started with seed i. We then evolve the solutions for 
another t — lOr-g. Because b32 has reached a saturated 
state long before the 85-th output, the energy levels re- 
main almost constant in all these simulations. The result 
is an ensemble of 16 saturated solutions [15]. 

The positions of the condensate vortices are simply 
given by the phases of the k\ modes. We unfold the 
trajectories from the computational domain [0, 2tt) 2 into 
M 2 and shift their starting points to the origin. We label 
respectively the x- and y-displacements, respectively, by 
Ax and Ay and plot such a trajectory (for path 100 ) in 
Fig. 5. The small gray box near the origin shows the 
size of the original domain [0, 27r) 2 . With the unfolded 
trajectories, we compute the displacement square Ar 2 = 
Ax 2 +Ay 2 for all 16 path 1 runs. For each simulation, we 
plot a gray curve in Fig. 6. The thick solid curve shows 
their (ensemble) average. In the following, we describe a 
way to model the motion of this condensate. 

We define the "mean" vorticity by 

ZJ(x)= ^e ikx (21) 

|fe|=fci 

so the fluctuating vorticity is u' = u> — ZU. Averaging the 
Navier-Stokes equation, we obtain 

d t uJ+ V • (uU) = -V -!F+ vV 2 U. (22) 

In the above equation, J 1 = ulu — uuJ is the space- 
dependent mean vorticity flux. Note that the forcing g 
is at small wave numbers fcj » k\, so its mean vanishes. 
It is straightforward to show that the non-linear term 
V • (uuj) vanishes identically. Both the creation and the 
movement of the condensate, therefore, must be due to 
the flux ~T. 

We now model this mean vorticity flux using the usual 
technique of mean-field theory 

Ti = aiZJ + PijdjUJ. (23) 

The transport coefficients Qj is usually referred to as the 
kinetic anisotropic alpha effect [16-18], and j3ij is a neg- 
ative eddy diffusivity tensor. Both at and are con- 
stants in space because we consider only the k\ modes. 
The mean-field equation then becomes 

(d t + aidi)oJ = (v6ij - /3ij)didjoJ. (24) 

The coefficient a, cannot change the amplitude of the 
condensate. It can be thought of as the Lagrangian veloc- 
ity of the pair of condensate vortices. The inverse energy 
cascade, therefore, must be described by anti- diffusion. 

Expanding V • JF, there are only two independent 
Fourier coefficients, T x m& an d -FyMyi * na * en ter the 
mean-field equation (the first subscript denotes compo- 
nent, the second one denotes wave vector). Comparing 
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FIG. 5: Unfolded trajectory of the positive condensate vortex 
in simulation path 100 . The computational domain is [0, 2tt) 2 
(see gray box near the origin). We first unfold the trajectory 
into R 2 and then shift the initial position to the origin. 



FIG. 6: The gray curves are the displacement square Ar 2 for 
all 16 path 1 runs. The thick solid curve shows their average. 
The dashed curve, which is indistinguishable from the solid 
curve, is the solution (35) with parameter £ = 0.1. 



to our parametrization (23), 



Re 
Im 



x,kix 
i 

J~ x.kix 



0. 



uu 



Re 



Im 



J~y,kiy _ 
J~y,kiy 



(25) 
(26) 



It is clear that (/3 XX ) — (/3 yy ) = v at saturation. 

The value of a, however, is much more difficult to ob- 
tain because there is no constraint by conservation laws. 
We can only perform a rough estimate: The flux T is a 
convolution in Fourier space. Given that kfE^/Z'^ > 1, 
the condensate-fluctuation interaction dominates, so 



u 



(27) 



where u^ ki represents the typical value of a band-pass 
filtered velocity. Comparing the estimate with Eq. (25) 
and employing the results from our three-scale model, we 
obtain 



(a 2 ) ~ fci£y5 fei 



fi 



2vk 2 



2/3 



4/3 



2kl /3 ' 



(28) 



2 (a 2 ) 



Note that the same argument leads to (/3 2 ) 
Fortunately, it does not contradict (0) = v. 

Let us now take a to be statistically independent of 
tJ. Equation (24) then becomes the equation of a passive 
scalar advected by a random velocity a. As a has a 
stationary variance (28), its simplest model would be an 
Ornstein-Uhlenbeck process [19-21] 



dtot = —a/i 



0. 



(29) 



In the above phenomenological equation, 1/t q is an ef- 
fective drag coefficient and <p i s an effective stochastic 



forcing. Because a is at a scale close to the condensate, 
we model the effective forcing <fi by Gaussian white noise 
with zero mean and variance equal to £e. Here £ is a 
tunable parameter and e is the mean energy input to the 
condensate. Standard Ito calculus gives 

(a(t)) = a(0)e-* /r °, (30) 

(a(s) • a(t)> = a(0) 2 e-( t+s )/ r = 

■„1 _ f,2 



uf ki 



Jkl-kf 



-(t-s)/r a 



-(t+s)/r Q 



(31) 



where the second equation is only valid for s < t. 



By requiring (a(i) 2 ) = (a 2 ) for arbitrary large t, we 



can solve for the time scale 

1 



7*2 



k\ 



2£/ j 2/3 fc 2/3 l k l 



(32) 



Furthermore, starting with an expected initial condition 
a(0) 2 = (a 2 ) allows us to simplify the correlation to 



(a(s) ■ a(t)) = (a 2 )e 



-{t-s)/r a 



(33) 



The displacement of the condensate can now be solved 
by using the simple equation 



d t r 



(34) 



The motion of an inertial particle in a fluid with Stokes 
drag is described by the same pair of equations [Eqs. (29) 
and (34)]. The role of effective velocity is played by a 
and the effective Stokes time is given by r a . 

Let Ar(i) = r(t) — r(0) be the displacement from the 
initial position, we have 



(35) 



(Ar(t) 2 ) = 2 / dt' / ds' (a(s') ■ a(t')) 
Jo Jo 

= 2(a 2 )r Q [i + T Q (e-'/ r ° - 1) 
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FIG. 7: We apply coordinate transformations to shift the 
center of the positive condensate vortex to the origin in sim- 
ulation c32. The noisy gray curves show the vorticity at sev- 
eral different snapshots along the a;-axis (with y = 0). The 
black solid curve is an average over 1000 such profiles. It con- 
sists of a sharp Gaussian core (dashed curve) and exponential 
wings (dotted curve), although for \x\ > tt/2 it falls off super- 
exponentially. The open circles show the reconstruction from 
the Fourier coefficients listed in Table II. 
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FIG. 8: The solid and dotted histogram are the energy spectra 
Ek of the full (oj) and the coherent (£j) vorticity in simula- 
tion c32 at t = 10000. The spectral bins have even width 
in log-scale. The valley and second peak right next to the 
condensate mode are physical. The valley corresponds to the 
fe = ^/2kl and 2fci modes, which have even k x + k y (open cir- 
cles in inset). The second peak corresponds to the |fc| = V5fci 
modes, which are the first harmonics of the condensate (gray 
circles in inset). 



Its asymptotic behavior is 



(Ar(i) 5 



(a 2 ) t 2 , t< T a : 
2(a 2 )r a t, t > r a . 



(36) 



The early time behavior is completely determined by 
(a 2 ), while the free parameter £ controls only the Brow- 
nian motion of the condensate at late time. 

The dashed curve in Fig. 6 is Eq. (35) with a parameter 
£ = 0.1, which corresponds to r Q — 6.384. It is practi- 
cally indistinguishable from the thick curve, except at 
very late times (At > 3000). Note that at short times 
the condensate movement is independent of £. Hence 
the agreement between the thick solid and dashed curves 
for At < r a provides support to our phenomenological 
model. Although we are not able to derive £, the late 
time agreement suggests that the condensate motion is 
Brownian for At > r Q . 



C. Shape of Coherent Vortices 

We can invert the procedure of computing the vortex 
trajectory to hold fixed the positive condensate vortex 
at the origin of the computational domain. In Fig. 7, 
the noisy gray curves show ten different profiles of the 
condensate along the x axis (with y = 0). Performing 
an average over 1000 such profiles, we obtain the black 
solid curve. To avoid confusion, we refer to this mean 
profile as the pair of coherent vortices and denote it by 
oj, while the name condensate vortices is kept for ul [see 



definition (21)]. The function oj(x, 0) can be well approx- 
imated by a sharp Gaussian core (fitted by the dashed 
curve) and exponential wings (fitted by the dotted curve) 
as shown in the figure. There is no preference between 
positive and negative vorticity so uj must change sign in 
the domain, which implies faster than exponential fall 
off for | a; | ~ L/2. This non-trivial shape cannot be fully 
described by the 4 modes with |fc| = ki. 

Nevertheless, it is possible to describe Q by a small 
number of Fourier coefficients. We find the Fourier coeffi- 
cients of the coherent vortices satisfying \ujk\ > 10 _3 |o}fc 1 1 
and list the "upper triangular" part in Table II. The 
shown values are real because of Hermitian symmetry. 
The lower triangle and other quadrants are given by 
tik^^y = &±ky ,±kv because of parity and discrete rota- 
tional symmetry. The open circles in Fig. 7 show the 
reconstructed profile from these coefficients. Despite the 
fact that we only have a few modes, the approximation is 
extremely good in both the exponential wings and the far 
tails where the power falls off more rapidly. The sharp 
Gaussian core, however, requires higher wave numbers to 
be represented correctly because of its short length scale. 

The omitted values in the table exhibit an interesting 
pattern. The Fourier coefficients with even k x + k y are 
significantly smaller than the ones with odd k x + k y . This 
is because the coherent vortices follow the same symme- 
try properties of the condensate, 

flip sign along diagonals 



j(x, y) = uj(x + L,y + L) = -Q(x + L/2, y + L/2). 



strictly periodic 



(37) 
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1.65 




0.79 




0.45 
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0.61 




0.32 
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2.54 




1.43 




0.72 




0.41 
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1.75 




0.8 




0.52 




0.14 
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0.6 




0.32 
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0.38 
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TABLE II: The "upper triangular" part of the Fourier coefficients of the coherent vortices, which are real because of Hermitian 
symmetry. The lower triangle and other quadrants are given by uJk x ,k y = &±k y ,±k x because of parity and discrete rotational 
symmetry. Coefficients that do not satisfy |a)fc| > lO^lu^ | are omitted in the table. We can reconstruct the shape of both the 
exponential wings and fast fall off in Fig. 7 (see open circles in the figure) by using the shown coefficients. The sharp Gaussian 
core, however, requires higher wave numbers to be represented correctly. The omitted values exhibit an interesting pattern 
that can be derived by symmetry, see Sec. Ill C for details. 



Fourier transforming the above equation, we immediately 
obtain the condition 

w fc .,fc„ =-£>*.,*, (-l) fc " +fc ». (38) 

Hence, the coherent vortices can only occupy modes with 
odd k x + k y . 

The above property has interesting implication for the 
energy spectrum of 2D turbulence. We graphically rep- 
resent the Fourier amplitudes in the inset of Fig. 8. The 
most energetic ki modes, which correspond to the value 
28.38 in Table II, are marked with filled circles. Their 
"higher harmonics", with value 10.37, are in dark gray, 
which form the second peak with radius y/Eki in the inte- 
grated spectrum as shown by the solid histogram. Other 
omitted modes with even k x + k y are marked as open 
circles, which form a spectral valley in the spectrum. 

The spectral valley and the second peak next to the 
condensate in Fig. 8 are, to our knowledge, not seen in 
previous numerical studies such as [5-7, 11]. It is possible 
that these earlier studies have not reached the late satu- 
rated stage we study, or perhaps because evenly spacing 
bins "wash" out the spectra as in Fig. 3. In the later 
case, reducing the bin size to ki/2 should make the spec- 
tral valley and second peak visible, although the spec- 
trum may become noisier. The spectrum obtained from 
atmospheric data [22] or from direct numerical simula- 
tions with Ekman friction [23] would also not have these 
features as in those cases energy from the large scales 
is removed by friction. Nevertheless, we should remark 
that in Fig. 3 of Hossain et al. [4], there is an indica- 
tion of such a second peak, although the authors did not 
comment on it. 

As we show by the symmetry argument, the non-trivial 
shape of the coherent vortices is responsible for these ad- 
ditional spectral features. To verify this, we plot the en- 
ergy spectrum of w by a dashed histogram in Fig. 8. The 
tight agreement between the solid and dotted histograms 
in small wave number strongly support our claim. 



IV. DISCUSSION 

In this paper, we use a three-scale model to derive the 
saturation time scales and saturation levels of the con- 
densate and turbulent fluctuation in 2D hydrodynamic 
turbulence. This requires integrating the 2D Navier- 
Stokes equations for a very long time. This has been 
made possible by virtue of the high performance of the 
CPUs. We use the Ornstein-Uhlenbeck process as a phe- 
nomenological model to describe the movement of con- 
densate vortices. In terms of the saturation time scales 
r-g- and Tz' , the DNS agree quite well with the analytical 
predictions. The saturation levels are, however, off by 
a constant factor £ = 0.76. We speculate that this dis- 
agreement is because our three-scale model ignores the 
non-trivial shape of the coherent vortices. The higher 
harmonics of the condensate modify the spectra at small 
wave number and affect the saturation level. We have 
checked that it is not possible to remedy this problem by 
just changing k\ in our model to an effective wave number 
&cond ^ k\. A more sophisticated model is needed. 

Before we conclude, we should point out that the direct 
(forward) enstrophy flux, 

%ir = Voo « fi(kf - k\ ), (39) 

is smaller than the enstrophy input f?k? in our three- 
scale model. There is no inconsistency here. Their differ- 
ence is simply the (very small) inverse enstrophy leakage 

r, inv = tfk? - Vdil « f? k{ (40) 

that we commented on for Eq. (9). The property ?7i nv <C 
?7dir is nothing magical. It is simply a consequence of the 
broken power law in our three-scale model. 

Similarly, we can check how the direct (forward) en- 
ergy leakage e^ir compares with the inverse cascade e; nv . 
Recalling that e is only the inverse energy transfer rate 
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at fci, we need to add the contribution from 2vZ' ', 



£di: 



4*7 In ^ 



# 1 



^4 

3 \ fc . 2 



(41) 



Using the original definition (20), it is easy to verify that 
1 — 37/2 < 1, which allows us to conclude einv ~> £dir- 
Most of the input energy is inversely cascaded and dissi- 
pated at the condensate. Note that the above approx- 
imations are made at the same order. The fact that 
£inv + £dir = S\ holds, therefore, shows the consistence 
of our model. 

Finally, our work stresses the fact that, as in three- 
dimensional turbulence [24] , the order of taking the limit 
of long time (t — > 00) and small viscosity [y — > 0) is also 
crucial in 2D turbulence. Specifically, taking the limit 
v — > first leads to a (linear) divergence in time 



lim 

t— fOO 



f lim e) = 



lim ftt. 

t—too 



(43) 



We must take the limits in the correct order, which is, 

f? Ikl 



lim ( lim E) 



lim 



lim 



2vk\ jkj - k\ 

8 



2vk{ 



(44) 



The system can then reach a steady state as long as 
the viscosity remains non-zero, even without the Ekman 
term. 
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